/*
 * Created on 08/08/2005
 *
 * TODO To change the template for this generated file go to
 * Window - Preferences - Java - Code Style - Code Templates
 */
package br.ufmg.dees.insane.solution;

import br.ufmg.dees.insane.util.IVector;

/**
 * @author jamile
 *
 * TODO To change the template for this generated type comment go to
 * Window - Preferences - Java - Code Style - Code Templates
 */
public class OrthogonalResidueStandardNewtonRapson extends StandardNewtonRapson {

	private static final long serialVersionUID = 1L;
	private IVector incXinicial;
    private IVector incX;
    private IVector incR;
    private double csi;
	private double dXmax;
	private double dIncFactor;
    
    public OrthogonalResidueStandardNewtonRapson() {
        super();
    }
	
	/**
	 * @param ass
	 * @param numMaxIterations
	 * @param tolerance
	 * @param convType
	 */
	public OrthogonalResidueStandardNewtonRapson(int numMaxIterations, double tolerance, byte convType, double dXmax) {
		super(numMaxIterations, tolerance, convType);
		this.dXmax = dXmax;
	}
	
	 public void assignStepState(double lambda){
	 	this.dIncFactor = lambda;
	 	incX.scale(this.dIncFactor, dxp);
	 	this.correct(incX);
	 	double aux1 = incX.dot(incXinicial);
        if(currentIteration == 1 && this.getLabel()==1){
            loadFactorTotal+=dIncFactor;
        }
	 	if(aux1 < 0.0){
	 		this.dIncFactor = -1.0 * this.dIncFactor;
	 		incX.scale(-1.0);
	 	}
	 	r.scale(loadFactorTotal,p);
	 	incR.scale(this.dIncFactor, p);
	 	x.add(x, incX);    
    	assembler.setXu(x);	
    }
	 
	 public IVector getUnbalancedX(){
    	f=assembler.getFp();
    	IVector aux1 = new IVector (numEq);
    	aux1.sub(r,f);
    	double aux2 = aux1.dot(incX);
    	double aux3 = incR.dot(incX);
    	csi = - aux2/aux3;
    	incR.scale(csi);
    	unbalancedForce.add(aux1,incR);
    	dx = sel.solverX(assembler.getIncrementalCuu(),unbalancedForce);
    	this.correct(dx);
     	return dx;
    }
	 
	 private void correct(IVector x){
    	double normX = x.norm();
    	if (normX > dXmax){
            x.scale(dXmax / normX);
        }
    }
	 
	 public void setDXmax(double dx){
    	this.dXmax = dx;
    }
	 
	 public double getCsi( ){
    	return this.csi;
    }
	 
	 public double getdIncFactor( ){
    	return this.dIncFactor;
    }
	
	public void execute(){    	
    	this.numEq = assembler.getSizeOfXu();
        this.incX = new IVector(numEq);
        this.incR = new IVector(numEq);
        this.incXinicial = new IVector(numEq);
		this.dxp = new IVector(numEq);
		this.dx = new IVector(numEq);
        this.r = new IVector(numEq);
        this.p = new IVector(numEq);
        this.f= new IVector(numEq);
        this.unbalancedForce= new IVector(numEq);
        double l_factor;
        this.x = this.assembler.getXu();
        
       	this.p = this.getProportionalLoad();
        this.r.scale(this.loadFactorTotal,p);
        this.r.add(this.getConstantLoad());
        
       	dx.zero();
       	incX.zero();
       	incR.zero();
       	incXinicial.zero();
    	currentIteration = 1;    			
    	for ( int j = 0; j < numMaxIterations; j++ ){
    		this.getXP();
            if ( currentIteration == 1 ){
                this.assignStepState(is.getPredictor());
                this.getUnbalancedX(); 
            }
    		if ( currentIteration == 1 ){
    			this.assignStepState(is.getPredictor());
    			this.getUnbalancedX(); 
    		}
    		else getUnbalancedX(); 
    		l_factor = is.getCorrector();
     		x.add(x, dx);
     		incX.add(incX, dx);
     		assembler.setXu(x);
    		setConvergence();
            if (!convResult) currentIteration++;
            else{
    			loadFactorTotal += l_factor;
    			incXinicial = incX;
    			return;
    		}
    	}
    }

}
